Lifetime of the first and second collective excitations in metallic nanoparticles 
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We determine the lifetime of the surface plasmon in metallic nanoparticles under various condi- 
tions, concentrating on the Landau damping, which is the dominant mechanism for intermediate-size 
particles. Besides the main contribution to the lifetime, which smoothly increases with the size of 
the particle, our semiclassical evaluation yields an additional oscillating component. For the case 
of noble metal particles embedded in a dielectric medium, it is crucial to consider the details of 
the electronic confinement; we show that in this case the lifetime is determined by the shape of 
the self-consistent potential near the surface. Strong enough perturbations may lead to the second 
collective excitation of the electronic system. We study its lifetime, which is limited by two decay 
channels; Landau damping and ionization. We determine the size dependence of both contributions 
and show that the second collective excitation remains as a well-defined resonance. 



PACS numbers: 78.67.Bf, 73.20.Mf, 71.45.Gm, 31.15.Gy 



I. INTRODUCTION 

The surface plasmon (SP) resonance is a very impor- 
tant collective excitation in metallic clusters»iii2i It is the 
dipolar vibration of the electronic center of mass with re- 
spect to the positive ionic charge, analogous to the giant 
resonance of nucleiji^ Since an external electromagnetic 
dipole field couples directly to the electronic center of 
mass, the photoabsorption spectrum of a metallic cluster 
is dominated by the SP. The lifetime of this collective ex- 
citation is a determining factor in the relaxation process 
studied in femtosecond pump-probe experiments.* 

The classical electromagnetic theory for a charged 
metallic sphere in the vacuum yields the energy Hlom of 
the resonance, with the Mie frequency ujm = ^p/\/3, 
where ujp — (47me^/me)^/^ is the bulk plasma frequency 
of the metal, n, e, and TOo being the electron density, 
charge, and mass, respectively. If the clusters are em- 
bedded in a matrix (of dielectric constant em) and/or we 
consider noble metal clusters (where the effect of the d 
electrons can be modeled by a dielectric function Cd) the 
Mie frequency takes the form ujm — (^p/V^d + Sem. The 
Mie frequency is close to the experimentally measured 
resonances. Such an agreement is not surprising since 
we deal with a collective excitation with a clear classical 
counterpart. Small red- and blueshifts with respect to 
wm have been experimentally observed in different phy- 
sical conditions, and various microscopic approaches have 
been developed to account for the frequency shifts, The 
most successful among them are based on a jellium de- 
scription (where the ionic positive charges are taken as 
a uniform background) and linear-response theory in the 
framework of the time-dependent local density approxi- 
mation (TDLDA):-5 

While it is difficult to measure the SP lifetime, nu- 
merous data for the linewidths of the absorption peak of 
ensembles of nanoparticles are available ji^SiSii but their 



theoretical analysis has proven to be quite involved. In 
principle, inhomogcncous effects arising from the disper- 
sion among the probed ensemble of clusters have to be 
separated from the properties of single particles. Intrin- 
sic effects depending on the bulk properties of the metal 
have to be separated from size-dependent properties of 
the cluster, and from the effect of the interaction with 
the local environment (matrix). In addition, the decay 
of the SP may follow different channels depending on the 
size of the cluster We calculate the Landau damping 
contribution to the linewidth (i.e. decay into particle- 
hole pairs), which dominates in the case of small clusters 
of radius a in the range 0.5-2.5 nm.-- For larger particles 
the Landau damping competes with radiation damping. 

Recent measurements of single-cluster optical absorp- 
tion have rendered accessible the optical properties 
of individual nano-objectsi^i2ii£ Most of the individual 
nanoparticles studied so far (in statio^ or dynamio^ se- 
tups) are too large to be in the Landau regime. However 
the linewidth of a single 2.5 nm radius gold nanoparticle 
has been determined lately>iS The possibility of overcom- 
ing the inhomogeneous broadening, and the application 
as biological markers, resulted in a renewed interest for 
the optical response of metallic clusters. 

Kawabata and Kubo studied the Landau damping of 
the SPfii and using linear response theory they proposed 
a total linewidth 

rt(a) = ri + r(a), 

with Pi a constant intrinsic width and T(a) inversely 
proportional to the particle size a. Barma and 
Subrahmanyamji^ and Yannouleas and Brogliai^ im- 
proved this calculation and proposed corrections to the 
behavior of P outside the regime /iwm/^f ^ 1, where e-p 
is the Fermi energy. They obtained 



r(a) = |^5(a 
2 fcpa 



(1) 
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where fcp is the Fermi wave-vector, and g a function of 
the ratio ^ = Hlom/^f- Numerical calculations within 
the TDLDA on free alkaline clustersi^ agree with this 
analytical result for 1.5 ^ a ^ 2.5 nm. For smaller radii, 
F shows a nonmonotonous size dependence superposed to 
the overall behavior of Eq. These shell effects arise 
from the electron-hole density-density correlations in the 
angular-momentum restricted density of states, and a 
semiclassical evaluation of F is in good agreement with 
TDLDA calculations. 

The experimentally observed nonmonotonous size de- 
pendence of the plasmon linewidth for charged alka- 
line metal nanoparticles in vacuum^ is consistent with 
the theoretical calculations. However, our calculated 
linewidths yield lower bounds for the experimentally 
measured ones, and the corresponding lifetimes are upper 
bounds of those found in real systems. Further measure- 
ments with smaller radii and a narrower size-distribution 
seem necessary to clearly establish the connection with 
the theory. 

Noble metal clusters embedded in inert matrices^ also 
exhibit a nonmonotonous linewidth for small a. However, 
a direct application of Eq. overestimates the smooth 
part of F.l^ This discrepancy motivates us to develop 
a refined theoretical description of the SP lifetime for 
the case of clusters with internal dielectric constant ed 
and surrounded by a dielectric medium with constant e^. 
The presence of an inhomogeneous dielectric environment 
leads to the modification of Eq. ^ . 

Under a weak initial optical excitation, only the first 
surface plasmon (that we simply denote "surface plas- 
mon" when there is no possibility of confusion) is ex- 
cited. With sufficiently strong initial excitations, we can 
also reach the second quantum level of the center-of-mass 
motion, known as the second (or double) plasmon. Such 
a resonance will be experimentally relevant provided its 
lifetime is sufficiently large (in the scale of cjj^^). The life- 
time is given by the anharmonicities of the center-of-mass 
system and by its interactions with the other degrees of 
freedom. Like in the previous discussion, the Landau 
damping is an important channel for the decay of the se- 
cond plasmon, but a new channel appears when 2hL0M is 
larger than the ionization energy: the ionization in which 
the cluster looses an electron into the continuumii^ Such 
a process was discussed in order to interpret the ioniza- 
tion of charged Najj clusters observed by Schlipper and 
collaborators.^^ 

We calculate the decay rates associated with different 
channels for the single- and double-plasmon states using 
a semiclassical approach within a mean-field description 
of the nanoparticle. Whenever it is possible, we verify the 
semiclassical approach by comparing to numerical calcu- 
lations. We characterize the size-dependent oscillations 
of the first and second plasmon linewidths for the case of 
free alkaline metals. In addition, we analyze the theoreti- 
cal difficulties in extending these calculations to the case 
of embedded and/or noble metal clusters and propose a 
way to overcome them. We also apply our semiclassical 



approach to the calculation of the ionization rate via the 
double-plasmon channel, and obtain results comparable 
with the experiments^ 

The paper is organized as follows: In Sec. we in- 
troduce the basic formalism for the photoabsorption and 
the SP linewidth. In Sec. lllll we present the semiclassical 
calculation of the single-plasmon linewidth, testing some 
of the approximations that we will use in the sequel. In 
Sec. IIVI we study the case of noble metal nanoparticles 
embedded in a dielectric medium and present the need 
to improve the existing theory for this case. In Sec. El we 
show a semiclassical description of the two main channels 
contributing to the decay of the double plasmon: Landau 
damping and ionization. Finally in Sec. I VII we draw the 
conclusions and the perspectives of our work. We relegate 
to the appendix a few technical, but important issues; in 
Appendix ^ we extend the standard calculation of the 
plasmon linewidth to the case where the cluster is made 
of a noble metal and/or is embedded in a nonreactive 
matrix. In Appendix^ we show how to take advantage 
of the spherical symmetry in semiclassical calculations 
like the ones of this paper, and how to recover some well- 
known results. In Appendix O we present the frequency 
dependence of the different plasmon linewidths. 



II. PHOTOABSORPTION AND PLASMON 
LINEWIDTH 

When the cluster is placed in an electromagnetic field 
with a wavelength much larger than its size^i^ the pho- 
toabsorption cross section is obtained from the applica- 
tion of the dipole operator on the ground state of the 
system: 



a(c.) = ^^|(/|.|0)r<5(na 
/ 



3c 



Ef), (2) 



where c is the velocity of light; |/) and Ef are, respec- 
tively, the many-body excited states and eigenenergies 
of the electronic system. The ground state is noted as 
|0) and its energy is taken as zero. In Eq. 10) , the pho- 
ton degrees of freedom have already been integrated out. 
In order to describe the electronic system, we consider 
a closed shell nanoparticle (perfectly spherical with a 
"magic number" of atoms) within a jellium model. The 
Hamiltonian representing N valence electrons in a uni- 
formly positively charged sphere of charge -l-A'^e is given 
by 



N r 9 
P 



2^ + 



1 ^ 



2 ir,, - r,- 



(3) 



with the single-particle confining potential 



C/(r) 



27rne a' 



3r' 



-Anne'' 



r < a, 
r > a, 



(4) 
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where n = N/V is the electronic density and V = ina^ /3 
the volume of the particle. The potential U is harmonic 
inside the particle and Coulomb-like outside it. 

The experimentally obtained photoabsorption cross 
section is dominated by the surface plasmon (SP) reso- 
nance at the frequency lom- The width T of this resonance 
can, in principle, be calculated from the eigenstates of 
H appearing in Eq. Q. However, this procedure is in 
general exceedingly difficult, and thus various approxi- 
mation schemes have been proposed^ii^ Among them, 
the TDLDA (time-dependent local density approxima- 
tion) is a numerical approach based on the local density 
approximation!^ We will use this numerical scheme as a 
check of analytical approaches that instructs us on the 
physical underlying mechanisms to the plasmon decay. 

Since we are in the long-wavelength limit where the 
field couples only to the electronic center of mass, a par- 
ticularly useful decomposition of the Hamiltonian ^ is 



H = Hr 



Hid + He 



Introducing the canonical coordinates (R, P), the (har- 
monic) centcr-of-mass Hamiltonian is given by 



Hr 



p2 



1 



2Nmc 



ifrei is the Hamiltonian of the relative coordinates and 
He expresses the coupling between the two subsystems. 
Introducing the standard position, momentum and low- 
ering operators 



NrricLUM 



2h 



Q 



we can write 



\/2NmchujM 



PQ: Q 



Hem — fu^M ^ I 0,qO-Q + 2 
Q=X,Y,Z 



It is difficult to handle H^ci and He in the general case. 
A notable exception is that of a confining potential which 
is not given by Eq. but which is harmonic for all r. 
We are then in the situation for which Kohn's theoremiS 
applies. It states that in a purely harmonic confinement 
potential, and with interactions only depending on the 
interparticle distance, the center of mass and the relative 
coordinates decouple (i.e., He = 0). The motion of the 
center of mass is that of a harmonic oscillator, with the 
characteristic frequency of the confining potential, inde- 
pendent of the electron-electron interaction. Due to the 
decoupling, the SP has an infinite lifetime. Kohn's the- 
orem gives us a first insight into the relaxation process 
of the SP: The Coulomb part of U (for r > a) leads to 
the coupling of the center of mass and the relative coor- 
dinates (i.e. He ^0), and translates into the decay of the 
SP. 

For the realistic situation i7c 7^ 0, it is useful to de- 
scribe iJrci and He within the mean-field approximation. 



where H,-ci can be expressed as 



Hmf — 



where Sa are the eigenenergies for the mean-field poten- 
tial V and cjj (cq) creates (annihilates) the one-body 
eigenstate |a). Consequently, the mean-field approxima- 
tion to He will be given by the change SV induced in the 
one-body potential V by a displacement Z of the center 
of mass. In Appendix ^ we show how to obtain SV in 
a self-consistent way from the electronic Coulomb inter- 
actions [Eq. (|A2|) ]. In second quantization, we can write 



He = 



hrrie 



2N 



az)cicf3, 



(5) 



where dap is the matrix element of the classical dipole 
field between two eigenstates of the unperturbed mean- 
field problem [Eq. 

The laser excitation induces an initial electronic state 
corresponding to a rigid shift (with a magnitude Z) of the 
unperturbed ground state. Within our separation for the 
degrees of freedom of the electronic system, such an ini- 
tial state can be written as a product of the ground state 
for the relative coordinate system and a coherent state for 
the center of mass (along the direction of the excitation). 
Since the amplitude of the perturbation is assumed to be 
small, the initial coherent state can be approximated by 
a linear superposition of the ground state |Ocm) and the 
first (harmonic oscillator) excited state |lcm,z)- The life- 
time of such an initial state is that of the SP. It is related 
to the transition rate F of |lcm,z) to |Ocm) by Ti = h/T, 
while the dephasing time is given by T2 = 2Ti. This de- 
cay is due to the coupling He and results in the transition 
of the relative coordinate system from its ground state to 
excited ones (that within our mean-field assumption we 
note |Omf) and |/mf), respectively). 

Assuming a weak coupling He, the SP linewidth can 
be obtained from the Fermi Golden Rule as 



r = 2^^|(0c 

/mf 



E/mf) 



According to form He, the final mean-field states 

I/mf) are particle-hole excitations, and therefore 



F = 



N 



ph 



(6) 



where \p) and \h) represent, respectively, particle and 
hole states of the mean-field problem. 

Form jni of the SP linewidth can also be derived from 
discrete-matrix random phase approximation^^ using the 
classical field associated with the collective state as the 
source of particle-hole transitions. The procedure pre- 
sented above is easy to generalize for the two cases impor- 
tant for our work: a nonhomogeneous dielectric function 
and the excitation of the second plasmon. 
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FIG. 1: Self-consistent potential as a function of the radial 
coordinate (in units of the Bohr radius oo) from the TDLDA 
calculations for a 832-atom nanoparticle with mean distance 
between electrons rs = (3/47rn)^''"^ = 3.03 ao, corresponding 
to a ~ 28.5 ao. The different curves are for e = ~ 
between 1 and 4, showing that the slope of the potential de- 
creases with increasing values of e. The corresponding Fermi 
levels are indicated by horizontal lines. 



terms of the Wigner-3j symbols as 



la 1/3 1 \ (la Ip ^ 



^ ' 6 j I j ■ 



The dipole matrix element of the radial problem can be 
written as 



7^' 



AV 



- m.is^-s,y Jo "^"'"^''^ dT 

(10) 

In the limit of large Vq, we have Uki{r) — 
^/2[a'^/'^ ji+i{ka)]^^r ji{kr) , where ji are the spherical 
Bessel functions and the allowed values of k are given 
by the quantization condition ji{ka) = 0, one obtainsi^ 



7^' 



(11) 



III. SIZE DEPENDENCE OF THE PLASMON 
LINEWIDTH 

In order to evaluate the plasmon linewidth from 
Eq. (jBJ, we need a description of the eigenstates |a) {p or 
h) of the mean-field problem. The self-consistent poten- 
tial obtained from TDLDA (thick line, Fig. [T} suggests 
that for analytical calculations, V{r) can be approxi- 
mated by a spherical well of radius a and finite height 



(7) 



with Q the Heaviside distribution. This stair like appro- 
ximation becomes more appropriate as the particle size 
increases. As we discuss in the next chapter, the dielec- 
tric constants inside and outside the cluster influence the 
steepness of V{r). 

The spherical symmetry of the problem allows us to 
separate the wave functions and matrix elements into ra- 
dial and angular components 



The summations appearing in Eq. ^ can be replaced 
by integrals provided one knows the particle (and hole) 
density of states (DOS). Decomposing the latter as a sum 
over its fixed angular momentum components [q{s) = 

Ei=o E™=-i we have 



47r?i 



Nnica'^UJM 7max(£F,?it^M) 



with Eh 



hujM ■ The overall factor of 2 accounts for 



the spin degeneracy. The angular part © of the dipole 
matrix element contains the selection rules ruh ~ nip and 
Ih — Zp± 1. Performing the sum over nip and Ih^ with the 



change of variables Sp 
and 77 = ka), we have 



F =4e^7 / dT]p r^l-qh'' 



(12) 



X ^ [ipQi^^ii'qh) + ilp + 1)^(^+1(77/1)] Ql^{Vp), 



and 



daf) — A 



Ida 



(8) 



Uki satisfies the radial Schrodinger equation IjBlj) . YJ™ 
represents the spherical harmonics, k = (2meey^^ /h is 
given by the principal quantum number, while I and m 
are the angular momentum quantum numbers. The an- 
gular part of the matrix element can be expressed in 



where 7 = 2TTh^ /SNmlcoMa^ , ??™ = t?f max (1, V?), 
^max ^ r]Fy^l + ^, ^ = hujM/£F and r]F = kpa. 

The SP linewidth depends on the Mixed DOS of the 
particles and holes. The asymptotic distributions of the 
zeros of the Bessel functions were used in Refs. and 10 
to obtain the leading behavior of F for the largest radii of 
the considered interval. Corrections, relevant for smaller 
radii, necessitate numerical or semiclassical approaches. 
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A. Semiclassical approach and smooth 
size-dependent component of the plasmon Hnewidth 

The semiclassical approximation to the radial problem 
(see Appendix IB|) allows us to write the Z-fixed DOS as 



Si{e) in 



(13) 



The classical action of the periodic orbit at energy e is 



, (1 + 1/2 

/ H — arccos 

2/ \ ka 



while its period is given by 



n{e) 



{kaf -{1 + l/2f 



eo{ka) 



and we note f the number of repetitions of the periodic 
orbit. Within the semiclassical approximation, the fi- 
nite height Vq of the self-consistent potential is irrelevant 
since the classical trajectories at a given energy are not 
sensitive to the shape of the potential above this energy. 

In the semiclassical approach of Ref.0, that we extend 
and improve in the sequel, it is natural to decompose the 
DOS into a smooth part and an oscillating part g"^"^ 
[Eqs. lO and (EU]. With Eq. lO, this leads to the 
dominant (smooth 1/a-dependent) component of F (due 
to the terms g'/ g^^ of the product) with nonmonotonous 
(in a) corrections. 

For the smooth part, we assume that Zp ^ 1, consis- 
tent with the fact that we are interested in leading order 
contributions in h. Then we use Z/j ± 1 ~ Zp and approxi- 
mate the sum over Ip by an integral. Setting y = Ip/rjp 
and z = ?7p/?7p, we find 



F°(a) = 



7(fcFa)' 
27r2 



dz 

max (1,^) 



(14) 

Performing the integrals of Eq. I|14|l leads to the smooth 
component F" given by Eq. The 1/a dependence 

agrees with the linear-response result of Kawabata and 
Kubo^ The function g appearing in Eq. |^ decreases 
with ,f with g{Q) — 1 and limf^op q(£) = 0. Its explicit 
form can be found in Refs.lT^ andll.ll: it is reproduced in 
Fig. of Appendix O 

The smooth component of the linewidth of the collec- 
tive state is inversely proportional to the radius of the 
nanoparticle: This has been interpreted^ as a surface 
effect arising from the confinement of the single-particle 
states. The analytical evaluation of F" agrees with the 
numerical calculations (see dashed line of Fig. |2J|. Ex- 
periments on charged alkaline clusters with a diameter 
in the range 1-5 nm in vacuum.® yield a linewidth of the 




FIG. 2: Inverse lifetime of the first collective excitation in Na 
nanoparticles as a function of the radius a of the particle. The 
dashed line is the smooth part of the single plasmon linewidth, 
Eq. The full line is the smooth part plus the oscillating 
contribution from Eq. II18II for a number of repetitions f — 1. 
This semiclassical result is compared to numerical TDLDA 
calculations (dots) for clusters with magic numbers of atoms 
between 20 and 1760. 



order of F ~ 1 eV. Although the charged character of 
those clusters limits the applicability of our model, we 
note that our calculated value is smaller, but of the same 
order of magnitude than the experimental one. This dif- 
ference might be explained by additional contributions 
to the linewidth present in the experiment. 



Shell effects and nonmonotonic behavior of the 
plasmon linewidth 



The oscillating part of the DOS H13|l gives rise to terms 
of the type gf g°^^ as well as g°^^ g°^^ . The former become 
negligible (in the semiclassical limit of small h) in Eq. I|12|) 
because one integrates a smooth function multiplied by 
a highly oscillating one. The latter yield 



Ip a—p,h 

' Si^iVa) 37ry 



E 



X > COS 



where //^ = Ip for 1^ ^ lp-1 and /;,^ = Ip + l for h = lp + 
1 . We can expand the product of the two cosines and keep 
only the contribution in leading order in ft, neglecting the 
highly oscillating term as a function of the particle and 
hole actions. We now write this contribution with the 
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help of the Poisson summation rule to obtain 



m— — oo fp,rh^l 
cr=± 



dL 



-1/2 



X E n + 1/2)2 e^'^ 

ih=ip±l a=p,h 

where we have defined the phase 



(15) 



h h 

37r _ „ ~ , 
- — (rp - r/,.) + 2TTmlp. 



(16) 



Performing a stationary phase approximation, given by 
the condition d4>ldlp\j — with the stationary points 

Zp, we obtain the stationary phase equation 



r„ arccos 



Ip + 1/2 



r/i arccos 



[ft + 1/2 



The phase of Eq. (|16|l indicates that the major contribu- 
tion to the integral over Ip in Eq. H15|l will be given by 
fp = fh and ifi — Q. We then select only one point within 
the full mesh of the stationary points, 



Tp + 1/2 _ + 1/2 



■Hp 



Vh 

1/2 



(17) 



Noticing that rjh — {rjp — VfO < Vpi we see that in 
order to satisfy Eq. (|17|l . we have to set Ih = Ip — 1. The 
stationary point is then given by 

J ^ Vp + Vh 
^ 2(77p - r]h) ' 

Performing the integral over Ip with the help of the sta- 
tionary phase approximation finally provides the follow- 
ing result for the oscillating part of the first plasmon 
linewidth: 



AT? 



dp' 



(3 + (3' 



: cos < 2f 



^{kFafifi-P'Y 



kFa{l3 - (3') 



(18) 



where /?' = yJP — S, and ^ — Itium/sf- The remaining 
integral over /3 can be performed numerically (solid line. 
Fig. 121). Assuming that fcpa ^ 1, using that fi — (3' ^ 
1 and that the sum over the number of repetitions is 
dominated by the first term, we see that the argument of 
the cosine is close to 2fcFa and 



cos (2/cFa)- 



Therefore the linewidth of the single SP excitation has 
a nonmonotonic behavior as a function of the size a of 
the metallic cluster. This is due to the density-density 
particle-hole correlation appearing in the Fermi Golden 
Rule ©. Let us mention that the result of Eq. H18|l is 
slightly different from the one of our previous work^ 
This is due to the fact that we have used here a more 
rigorous treatment of the semiclassical radial problem. 
As in Ref. 0i we have to set a phase shift in our analyti- 
cal prediction of Eq. 118|l to map the TDLDA numerical 
points in Fig. |21 This is due to the fact that we have 
taken only one stationary point [Eq. (|17|l ] and neglected 
all the other contributions coming from the full mesh of 
stationary points which influence the phase appearing in 
Eq. 113). 

A nonmonotonic behavior has also been observed ex- 
perimentally in the case of charged lithium clusters^ 
Our numerical TDLDA calculations (that we have ex- 
tended here to larger sizes, up to 1760 atoms) confirm 
the presence of size-dependent oscillations for alkaline 
metals. The semiclassical approach also predicts a non- 
monotonous behavior of P for noble metal clusters, in 
agreement with recent experimental resultsi^S However, 
the presence of different dielectric constants inside and 
outside the cluster render the problem more involved. 
This issue is discussed in the following section. 



IV. PLASMON LINEWIDTH WITH AN 
INHOMOGENEOUS DIELECTRIC 
ENVIRONMENT 

In a previous analysis of the surface plasmon (SP) 
linewidthfi^ we were interested by the case of noble 
metallic nanoparticles (where the d electrons are mod- 
eled via a dielectric constant ed) embedded in a matrix 
of dielectric constant em- The two dielectric constants 
affect lum as discussed in the introduction. However, a 
generalization of the derivation of Sec. IHII shows that, 
as long as we work with the hypothesis of a steep po- 
tential [Eq. (|7|)], the smooth part of F is still given by 
Eq. (Q. For a silver nanoparticle (ed — 3.7) embedded 
in an argon matrix (cm — 1.7)ri using Eq. (Q) yields a 
value of r° about three times larger than the TDLDA 
calculationsi^ (themselves in good agreement with exist- 
ing experiments'"). This discrepancy makes the more sys- 
tematic study of the dependence of the plasmon lifetime 
on ed and em presented in this section necessary. 

In Fig.|2Ia) we present the SP linewidth obtained from 
TDLDA for several particle sizes between N = 138 and 
1760, taking ed = 4 and e,ii — 2 and the electron density 
of silver (rg = 3.03 oq). As in the case of Fig. |21 we 
see that for relatively large radii the linewidth can be 
approximated by F" — C/{a/ao) while for smaller radii a, 
superimposed oscillations become noticeable. As shown 
in Fig.|2Ib), when plotting the coefficients C as a function 
of ed and em, we see that the numerical results are at odds 
with the simple prediction of Eq. ^ (upward continuous 
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FIG. 3: (a) Surface plasmon linewidth from the TDLDA as 
a function of the inverse radius for the example of td = 4 
and Eni = 2 (dots). The straight line is a linear fit F = 
C/(a/ao). (b) Prefactor C of the smooth 1/a size-dependent 
component of the surface plasmon linewidth F° as a function 
of Em for Ed = 1 (solid line), ed = 2 (dashed line), Ed = 3 
(dashed-dotted line), and ed = 4 (dotted line). The crosses 
connected by straight lines (guide-to-the-eye) represent the 
TDLDA calculations, while the increasing curves in the upper 
part of the figure depict the analytical expression Q. The 
thin gray line is for ed = em- The results presented in the 
figure are for the electron density of silver (rs = 3.03 ao). 



the self-consistent potential of a nanoparticle consisting 
of = 832 atoms (rg = 3.03 oq) for various values of 
e = Ed = Cm- This choice does not correspond to a physi- 
cal realization, but it is useful for the interpretation of the 
analytical work, as it merely represents a renormalization 
of the electronic charge. The main effect of increasing e is 
the decrease of the slope of the potential near the bound- 
ary r = a. This indicates that our approximation of a 
square-well potential becomes less valid as the dielectric 
constant is increased. The F*' dependence on e in this 
case is obtained by moving along the line ed = Cm hi 
Fig.Hb). 

In the following we refine the calculation of the dipole 
matrix element (jlOII in order to take into account the be- 
havior of the slope of the self-consistent potential. The 
finite value of the slope of the self-consistent potential is 
often ignored. But here, it is necessary to go beyond the 
hypothesis of infinitely steep potential walls [Eq. Q] in 
order to make progress. As it can be seen from Eq. UlUfl. 
the dipole matrix element is proportional to the matrix 
element of the derivative of the potential V with respect 
to r. In the sequel, we show that below a certain value, 
the dipole matrix element is directly proportional to the 
slope of the self-consistent potential near the interface 
and estimate the slope from a simple model. Since the 
linewidth is proportional to the square of the dipole ma- 
trix element, we see that F decreases with the slope, and 
thus with the increase of the dielectric constant. 



A. Surface plasmon linewidth with a soft 
self-consistent potential 

In order to improve our understanding of the role of 
a dielectric mismatch to the SP linewidth, we now need 
to come back to the evaluation of Eq. ^ without mak- 
ing the approximation of an infinitely steep well for the 
self-consistent potential. A simplified way of taking into 
account the noninfinite slope of V{r) is to change Eq. ^ 
by 



curves). 

The increase of F° with ed and em in the latter case 
arises from the fact that the function g is decreasing with 
^ = Hlom/sf and the Mie frequency ujm — Wp/V^d + ^e^n 
is redshifted when ed or em is increasing. Calculations 
performed for the electronic density of sodium (rg = 3.93 
ao) give the same kind of discrepancy between Eq. 
and TDLDA results-^S 

The discrepancy between the numerics and Eq. 
shows that a direct application of the analytical approach 
presented in the preceeding section does not reproduce 
the TDLDA results. As we will see in the following, the 
discrepancy is caused by approximating the electronic 
self-consistent potential by a square well. 

The TDLDA calculations show that the shape of the 
self-consistent potential is modified when one increases 
the dielectric constants ed or Cm- In Fig. ^ we present 



0, 



V{r) = < 



s r — a — 



r < a - — , 
2 

ds dg 
+ Vo, a - y r < a + y, 

ds 

r>a+-, 



where the distance d^ on which the slope s — Vo/dg 
is nonvanishing is assumed to be small as compared to 
the nanoparticle radius a. We first need an approxi- 
mation for the dipole matrix element between particle 
and hole states in that potential. As explained in Ap- 
pendix IB 31 this can be done semiclassically using the 
limit in which particle and hole states are close in en- 
ergy [{sp — eh)/eF — hajM/ep ^ !]• This semiclassical 
approximation relates the dipole matrix element to the 
Fourier components of the classical trajectory in the one- 
dimensional effective potential V°^(r). As a simplifying 



8 



approximation, we neglect the centrifugal part of the ef- 
fective potential above r > a — ds/2. Integrating the 
classical equation of motion, we obtain periodic trajecto- 
ries (for e < Vb) given by 



rit) 




Vn-e 



t S$ tr 



t > u 



with r{tc) — a — ds/2 and where t/ is the period. We 
can now evaluate the dipole matrix element using the 
semiclassical Eq. (|B7|I , neglecting the acceleration of the 
particle for r — a -I- ds/2 0~ (justified for a ^ dg). An 
expansion in 1/An (where An is the difference between 
the radial quantum number of the particle and of the 
hole) gives, up to an irrelevant phase factor 



ft3 



Tip [ep -EhY 



sin nAn 



(19) 



with 6t — Ti /2 — tc the time spent by the particle in the 
region where the slope is nonvanishing. 

An estimation of the argument of the sine gives 
{hiOM./ A){ds/a), with A the mean level spacing. Typ- 
ical values give Hujm/A ~ 10** » 1. In the limit of a very 
large slope, ds/a tends to zero. Then, the argument of 
the sine is very small compared to one, and we recover 
the semiclassical evaluation of Eq. I|B8I) with an infinite 
slope. On the contrary, if we assume that ds is of the 
order of the spillout lengthS (^ ao), the argument of the 
sine is much greater than one. Inserting Eq. (|19|l into 
Eq. lO, we obtain 



2s2 



Ae 



nAn- 



St 



Averaging the highly oscillating sine (squared) by 1/2 
gives for the SP linewidth in the limit ^ ^ 



r"(a) 



3s^ 



1 



1 



4 mnW?, fcpa 



(20) 



We then see that in the case of a soft self-consistent po- 
tential, the SP linewidth is proportional to the square 
of the slope s of that potential. When one increases the 
dielectric constant of the medium, the slope decreases 
(see Fig. and therefore r° decreases. We also notice 
that the smooth 1/a size dependence of the SP linewidth 
remains valid even for a finite slope. 



B. Steepness of the self-consistent potential with a 
dielectric mismatch 



In order to estimate the slope of the self-consistent 
potential, we consider the simpler geometry of a metallic 



slab of dielectric constant ed, bounded by two interfaces 
at a; = ±w/2 and with an infinite extension in the {y, z) 
plane, surrounded by a dielectric medium with a constant 
Em. This geometry allows us to simplify the problem to 
an effective one-dimensional system and can be expected 
to provide a good approximation for the shape of the 
potential near the interface for the sphere geometry. 

We make the jellium approximation for the ionic den- 
sity ni{x) = niQ{w/2 — |x|), with Q the step function, 
and work within the Thomas-Fermi approach, writing 
the local energy in the electrostatic field (j> as 



2mr 



and the electronic density (at zero temperature) as 



2mc 
If 



3/2 



3/2 



with fj, the chemical potential in the potential V{x) = 
—e(j){x). The Thomas-Fermi approach to surfaces is 
known to have serious shortcomings^^ (for instance, it 
predicts a vanishing work function). However, it will 
be useful for our estimation of the slope of the mean 
field seen by the charge carriers. The self-consistency is 
achieved through the Poisson equation 



d20 

dx^ 



■ 47re 
47re 



[nc{x) - Tii] , |a;| < 



nc(a;), 



|x|> 



(21) 



First we consider the simpler case where Cd = em = e- 
In this case, integrating once Eq. H21|l and invoking the 
continuity of the potential and the electrical field, we find 
for the slope of the self-consistent field at a; = w/2 



4e 



2mc 



3/4 5/4 
PI 
,7/4 



5e3/2 



3/2' 



-1 5/4 



where we have assumed the scaling /i w Mi/^ with 
the chemical potential in the case where e = 1, and e-p is 
the Fermi energy of the free electron gas. The chemical 
potential is fixed by the consistency condition 



du 



w 
2" 



(23) 



with 



If we do not have any dielectric constant (i.e., e — 1), 
the same equation is obtained but without the prefactor 
-^e. The integral in Eq. H23|l is clearly dominated by its 
prefactor. Then, assuming that the integral appearing 
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in this equation does not change appreciably when we 
have a dielectric constant, we find the scaling ^ « Mi/e- 
Therefore, we see that the slope at the interface is de- 
creasing with increasing values of the dielectric constant 
e, a feature confirmed by our TDLDA calculations (see 
Fig.[l|. 

In the case where we have a dielectric mismatch be- 
tween the metallic slab and the environment, the conti- 
nuity of the normal component of the displacement field 
D gives, perturbatively, in the limit jej — £^1 0, 



4e 





1 - 




f Ml 


ft2 ) 1/2 5/4 




^V^F 


/ \ 3/2 
Ed - Em / A^l \ 


1 - 


2 




Uf; 


K 3/2 





3/2 



3/2' 



1 5/4 



(24) 

with the scaling fi « //i/ed, which can be justified in 
the same manner as for the case of a single dielectric 
constant. The only difference is that in the case of a 
dielectric mismatch, we obtain Eq. (|23|l . up to a change 
of e by Ed- We then see that the slope s of the confining 
mean-field potential at the interface is decreasing either 
with Ed or Cm (for small |ed — eml), in agreement with the 
TDLDA calculations. 

This Thomas- Fermi approach to the mean-field poten- 
tial of a metallic slab then provides an estimate of the 
slope of that potential near the interface between the 
slab and the surrounding environment. It can be ex- 
pected that these results are also applicable to the more 
involved problem of the metallic sphere, up to some ge- 
ometrical prefactors. In the following section, we will 
incorporate our estimate of the self-consistent potential 
slope in our evaluation of the SP lifetime. 



Surface plasmon linewidth with a dielectric 
mismatch 



We can now use our estimate (|24ll for the slope of the 
self-consistent potential in our evaluation H2U|) of the SP 
linewidth. In order to do that, we assume that the chem- 
ical potential fii for e = 1 is the Fermi energy e-p of a free 
electron gas. 

In the case where we have a charge renormalization 
(i.e., ed — Cm = e), we obtain by inserting Eq. (|22ll into 
Eq. 1201) 



r°(«) 



9 ep 1 
5 kpa e^/2 



5e3/2 



5/2 



(25) 



This result qualitatively reproduces the decrease ob- 
tained from TDLDA for F" a/ao as a function of the di- 
electric constant e as it can be seen on Fig. [SJb) (gray 
thin line). We notice that for e = 1, we have F" « ep/kpa 
in the limit of small ^, which has to be compared with 
Eq. ^ giving 1.5 SF/k-pa. This small discrepancy is 



not surprising, regarding the various approximations we 
made here. 

In the case where we have a dielectric mismatch, by 
inserting Eq. H24(l into Eq. H2()|l and making the expansion 
for small Ae = ed ~ em, we obtain 



r°(a)=^rO^,=o)(«)+^Ae 



for fixed Cd and 



r°(a) 



BAe 



(26) 



(27) 



for fixed em. In the above two equations, A and B are 
two positive coefficients not specified here, and F^'^^^p^ 
is given by Eq. (|25|l . These results confirm the behavior 
of the TDLDA calculations depicted on Fig.|3Ib) around 
Ae = (thin gray line). For instance, if we are at em 
fixed, we see that when Ae > 0, Eq. (|27|l predicts that 
F° a/ao decreases for increasing value of ed. 

We have shown in this section how to take into ac- 
count an inhomogeneous dielectric environment in our 
semiclassical model through the corrections in the slope 
of the mean-field potential. This improved theory is in 
qualitative agreement with the TDLDA calculations. 



V. HIGHER COLLECTIVE EXCITATIONS: 
DOUBLE PLASMON 

In this section we discuss the lifetime of the second 
collective excitation level in metallic nanoparticles. Al- 
though there is no clear direct experimental observation 
of a double plasmon in metallic clusters, the development 
of femtosecond spectroscopy will certainly allow for de- 
tailed studies in the near future. Recent experiments 
observed the ionization of the charged cluster Na^jg by 
a femtosecond laser pulse and claimed it was a conse- 
quence of the excitation of the second plasmon state>ii'^ 
However, the analysis of the distribution of photoelec- 
trons yielded a thermal distribution and therefore the 
relevance of the double plasmon for this experiment is 
not yet settlediii'^ On the other hand, it is clear that a 
strong-enough laser pulse will excite the second collective 
state. Such an excitation will be a well-defined resonance 
only if its linewidth is small compared with other scales 
of the photoabsorption spectrum (like for instance lum)- 

Second collective excitations have been widely analy- 
zed in the context of giant dipolar resonances in nucleiiS^ 
The anharmonicities were found to be relatively small, 
making it possible to observe this resonance*?^ The the- 
oretical tools developed in nuclear physics have been 
adapted to the study of the double plasmon in metal- 
lic clustersi2LS& In particular, a variational approaches 
showed that the difference between the energy of the dou- 
ble plasmon and 2hu;M decreases as N~^^^ with the size 
of the nanoparticle. In our calculations, we will assume 
that the double-plasmon energy is exactly twice the Mie 
energy. 
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For most of the clusters of experimental interest, 
2Hlum > W > hujM, where W is the work function. Ioni- 
zation then appears as an additional decay channel of the 
double plasmon that competes with the Landau damp- 
ing, while it is not possible if only the single plasmon is 
excited. 



A. Second plasmon decay: Landau damping 

In this section, we consider processes which do not lead 
to ionization, that is, the final particle energies verify 
Ep < Vq = ep + W . A sufficiently strong laser excitation 
gives rise to an initial center-of-mass state which is a 
linear superposition of the ground-state |Ocm): the first 
(|lcm,z)), and the second (|2cm,z)) harmonic oscillator 
excited states. 

The second plasmon state can decay by two distinct 
Landau damping processes. A first-order process, with a 



rate r2^i, results from the transition of |2cm.z) (double 
plasmon) into |lcm,z) (single plasmon). The correspond- 
ing matrix element of the perturbation He between these 
two states is a factor of V2 larger than the one worked 
in Sec. ^ and then r2_,i = 2r [where F is the single- 
plasmon linewidth given by Eq. © and calculated under 
certain approximations in Sec. IIII| . Thus, the contribu- 
tion of the first-order process to the linewidth is just twice 
that of the single plasmon, and shows the same nonmono- 
tonic features superposed to a 1/a size-dependence. 

The other mechanism one has to take into account is 
the second-order process, where the double plasmon de- 
cays directly into the center-of-mass ground state. This 
is possible provided that Vq > 2hujM- In order to simplify 
the calculation we assume, for the remaining of this sec- 
tion, that Vo — > oo. The corresponding linewidth F2^o 
is given by the Fermi Golden Rule in second order in 
perturbation theory byi& 



2^0 = 2^^ 



(0cm, /MF|^?c|lcm,Z, /mf) (lcm,Z , /mF I -^c 1 2cm,Z , 



MF) 



M 



(5(2nwM-e/MF)- 



Expliciting the perturbation (j^l and restricting ourselves 
to the random phase approximation which allows only 
one particle-one hole transitions, wc obtain 

r2^o = ^^^^|i^p,,|'(5(2fi^M-ep + £/^), (28) 



Ar2 



ph 



with 



fu^M - £i + £h 



(29) 



The sum over i runs over all the virtual intermediate 
states. We use the same notations as in Sec. Illll and re- 
place the sums over particle and hole states by integrals 
over the energy with the appropriate density of states 
(DOS), which is approximated by its semiclassical coun- 
terpart. 

As in the case of the single plasmon, we work in the 
limit S> 1 in order to find the smooth size-dependent 
contribution F^^q (and the corresponding K^^y Using 
Eqs. (|H1), ©I Ullll - and the selection rules, we have 



K° - 



trripmh 



(30) 



Here, we have defined the integral 



drii 



Inserting Eq. H30|l into Eq. (|28|l . replacing the sum over 
Ip by an integral, we find 



^2^0 =2<r / drip rip 

J I'lP max ( 1 , ) 

dip Ip^ril - l^^rih"^ - II [lip{r]p,rih)Y , 



where the factor of 2 accounts for the spin degeneracy. 
We have introduced <; = QA{fROMY' / '5>t^'^ N'^ £q and rih = 
(?7p — 2r/p^)^/2. With the change of variables z = rip /rip, 



y — ^p/Vfj ^^"^ ^ — Vi/Vp' obtain 



2-tO 



(a) 



81 



10^3 (fc. 



■pa 



• HO, 



(31) 



where the function ft,(^) of the parameter ^ = /io^m/ef 
is smoothly increasing with h(0) = 0. An approximate 
expression of h is given in Appendix lO 

The total linewidth of the Landau damped second plas- 
mon is the sum of the first- and second-order processes: 
Fdp = F2_^i + F2^o- The different (smooth) size de- 
pendence of both processes [vp /a for the former and 
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{kFa)~^vp/a for the latter] implies that, except for the 
smallest clusters, the second-order process gives a negli- 
gible contribution to the linewidth of the double plasmon 
(in comparison with that of the first order). We might 
ask the question of whether the inclusion of the oscillat- 
ing components of both linewidths can affect the above 
conclusion in this range of particle sizes. An extension of 
the calculations presented in Sec. IIII 51 shows'^" that the 
oscillating part of the second-order channel of the double 
plasmon is given by 



Vir 



T^OSC / \ 



cos {2kpa). 



As indicated before, V^^j^i is given by twice the result 
of Eq. p8|l . therefore these nonmonotonic contributions 
cannot lead to a significant modification of our con- 
clusion about the irrelevance of the second-order term 
for the sizes of physical interest. We also notice that 
Tdp ^ since for typical nanoparticles, ep ~ huM 

and kpa 3> 1. Therefore, the Landau damping is not ca- 
pable of ruling out the second plasmon as a well-defined 
resonance. 

The lifetime of the second plasmon for the Landau 
damping processes is simply hT^p. From the experimen- 
tal point of view, what is usually more relevant is the 
time it takes for the double excited state of the center- 
of-mass system to return to its ground state rather than 
the lifetime of the excited state. Therefore, we also have 
to take into account the decay of the first plasmon into 
the ground-state Fi^o- If we assume that the recombi- 
nation of particle-hole pairs (created by the decay of the 
double plasmon into the single plasmon) is very fast as 
compared to other time scales, we have Fi^o — T- E)ue 
to the fact that lifetimes are additive, we have for this 
sequential decay a lifetime T2^i^q = 1.5hT^^. 



B. Second plasmon decay and ionization 

We now examine the last decay channel of the second 
plasmon state: the relaxation of this collective excitation 
by ejecting an electron from the nanoparticle (ionization, 
see inset of Fig. . We now need to determine the par- 
ticle and hole states in the self-consistent field [Eq. 
which has a finite height Vq, since the ionization process 
requires the states of the continuum. For simplicity, we 
will neglect the Coulomb tail seen for r > a by electrons 
with an energy Ep > Vq. 

In order to determine the particle and hole states, we 
close the system into a spherical box of radius L 3> a 
to quantize the states above the well and take the limit 
of L ^ oo at the end of our calculations. We need to 
do some approximations in order to simplify this diffi- 
cult problem. First, in the high energy limit, we assume 
that fcr 3> 1, and then use the asymptotic expansions 
of the quantum mechanical single-particle states inside 
and outside the well. Even though this approximation 
strongly affects the wave functions near r = 0, its impact 
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FIG. 4: Ionization linewidth of the second plasmon state 
as a function of the nanoparticle radius for singly charged 
Na clusters. Square: experimental value for Na^g taken from 
Ref.El We have assumed a constant work function W = 4.65 
eV and took the experimental Mie frequency of 2.75 eV. Inset: 
scheme of the ionization process of the double-plasmon state 
which decays by creating a particle-hole pair of energy 2hujM, 
via the intermediate state Since the energy of the particle 
is such that £„ > Vb, ionization occurs. 



on the dipole matrix elements is very small^ Second, 
for the states with energy e < Vq, we neglect the expo- 
nential decay of the wave function for r > a. Finally, 
in the spirit of the scattering theory, we use a simplifed 
expression for the normalization of the free states above 
the well. The above assumptions result in the following 
radial wave functions inside the well (e < Vq) 



sin (fcr — /7r/2), r ^ a, 



0, 



r > a. 



The wave- vector k = {2mceY^^ /h is given by the quanti- 
zation condition ka = Itt/2 + rnr, with n a non- negative 
integer. Outside of the well (e > Vq), we have 

> , , I ^ J ai{k)A{k) sin (fcr — ^7r/2), r ^ a, 
""kM - \l L \ sm[K{r - L)], r > a, 

with K = {k'^ — 2moVb/ft^)^/^. We have introduced the 
abbreviations 




ai{k) = sign 



sin [K{a — L)] 
sin {ka — It:/2) 



A{k) = Wsin-^ [K{a-L)] 



— I cos^ [K[a — L)]. 



The ionization rate of the double-plasmon state Fjon 
is given by Eq. H28|) in the case where the final particle 
states p of the sum are in the continuum. Since the ef- 
fective (second-order) matrix element Kph [Eq. is 
given by a sum over intermediate states i, we now have 
contributions from cases where i lies in the well as well 
as in the continuum. 
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When i represents a state in the well, using the angu- 
lar momentum selection rules, we can write in the limit 
fcr > 1 



and 



( 1 \'f^i — '^h 



a Q;;p(fcp)A(fcp) 



(32) 



■p2 



COS (Akpiu) 



sin (A/cpja) 



(33) 



where = ka — kfs (a, P = p,h,i). 

When i represents a state in the continuum, T^-jj'^!^ 
can be obtained by exchanging [p <-> i) and (i <-> /i) in 
Eq. H33|) . For the remaining case, we have 



Iph _ g Q;;p(fcp)Q:;^ + i(fci)A(fcp)A(fcj 



COS (Afcpio) 



sin (Afcpifl) 



(34) 



^i,,ip±i + B{kp,ki), 



where 
B{kp^ ki 



cos (AKpiL) — sin (AKpiX)ci(| A^pija) 
AKpia[cos (AKpjL)sign(AKpi)si(|AKpi|a)] L 



with si and ci the sine and cosine integral functions. 

The semiclassical Z-fixed smooth DOS can be approxi- 
mated by 





-1/2)2 


{kaY 


y/{KLy-{H 


h 1/2)2 



{nay 



e < Vo, 
e>Vo. 



There is an obvious divergency that occurs in the sum 
of Eq. (|29() for Si = Sh, as it can be seen on the matrix 
element H32(l . However, a careful analysis shows that the 
contribution around that divergency vanishes because of 
the alternating sign when one integrates over n^. For 
£i = £p, there is no divergency in Eq. 1)34(1 . Therefore the 
dominant contribution to Kph is given by the divergency 
of the term l/{hujM — + £h) that occurs for < Vq in 
the regime we are interested in {fkuu < W < 2?iwm)- We 
then have for the ionization rate 

i2 V S{2huju -£p + £h) 



P>Vq 

h<€F 



d/T) 7 d ■ 



f^M - St + Sh f^M - £7 + £h ' 



where the factor of 2 accounts for the two spin channels 
and the dap are given by Eq. ((SJ with the approxima- 
tions ((32|l and 1(33(1 for the radial matrix elements. Fur- 
thermore, we can distinguish in the above equation two 
contributions: off-diagonal terms (i ^ j) which have di- 
vergencies of the principal value type and that we neglect 
here, and diagonal terms (i — j) yielding divergencies 
which determine Fjon- We smooth out the energy Si ap- 
pearing in the denominator by introducing an imaginary 
part of the order of the mean level-spacing 



A 



3ne 



3/2 



y/hwM + £h 



at an energy hujyi -f Eh- This standard procedure of 
smoothing the divergencies is of critical importance, and 
that is why in Ref. 16 the final result is presented as a 
function of A. Summing over li and to^, the remaining 
sum over the radial quantum numbers Ui can be done 
with the help of 



E 



fvjj 



M 



4A huj 



M 



■ £h 



For the smooth terms of the sum, we have taken their 
values at the divergency to obtain 



120N^L 



E^. / 



de, 



Qip{sp)Qip{sh) 



P max (Vb:2?iajM) 

[A{kp)? 



{fuuJM + eh)WhujM +£h- ^/Sh) V^WM + £h) 



with Eh — £p — 2hwM- 
Taking the limit of L 



we finally obtain 



rion(a) ^ ^-n^<ii^X), 

80 Kpa 



(35) 



where ^ = hajM/^F and ( = W/sp. The function q of the 
two variables ^ and C is defined in Appendixinj Eq. l(Cip . 

The size scaling of Fjon is mainly given by a 1/a- 
dependence of the prefactor (Fig. 0J), despite the fact 
that the work function W appearing in the parameter 
C is size dependent and scales (for a neutral cluster) as^ 
W — Woo + 3e2/8a where Woo is the work function of 
the bulk material. 

Using the work function W — 4.65 eV and the exper- 
imental value of Hlom = 2.75 eV for the charged Na^Jg 
clusters of Ref. [13, Eq. ^33) yields Fjon ^ 0.1 eV, which 
corresponds to an ionization lifetime of the second plas- 
mon of 6.6 fs. This value is of the same order of magni- 
tude as the experimentally reported lifetime of 10 fs. It 
is also in rough agreement with the estimation yielded by 
the numerical calculations of Ref. 16^ based on a separa- 
ble residual interaction (10 to 20 fs). Therefore, despite 
the approximations we have been forced to make in our 
analytical calculations, we believe that we kept the es- 
sential ingredients of this complicated problem. The life- 
times obtained by the different procedures consistently 
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establish the second plasmon as a well-defined resonance 
in metallic clusters. While the numerical calculations of 
Ref. have been performed for just one size, our re- 
sults exhibit a clear size dependence that can be tested 
in future experiments. 

VI. CONCLUSION 

In this work we have analyzed the lifetime of collective 
excitations in metallic clusters. Different decay mecha- 
nisms have been studied within a semiclassical approach 
for the mean-field self-consistent potential describing the 
electrons in a jellium background. We have consid- 
ered Landau damping, which is the dominant relaxation 
mechanism for nanoparticles with radius a in the range 
0.5-2.5 nm. We found that the linewidth of the single sur- 
face plasmon exhibits a 1/a dependence, superimposed to 
an oscillating behavior arising from electron-hole density- 
density correlations. These results are in good agreement 
with numerical time-dependent local density approxima- 
tion calculations, and consistent with experiments on free 
alkaline nanoparticles. 

To describe noble metal clusters, we have taken into 
account the screening efi^ect of the d electrons and the 
modifications induced by the dielectric properties of an 
eventual matrix. We have demonstrated that such an in- 
homogeneous dielectric environment of the nanoparticles 
strongly affects the steepness of the self-consistent poten- 
tial, which in turn has a crucial influence on the plasmon 
linewidth. We could then solve the discrepancy presented 
in Ref. between the well-known Kawabata and Kubo 
formula on one side, against experiments and numerical 
calculations on the other side. The size-dependent os- 
cillations of the linewidth also depend on the dielectric 
constants through the slope of the self-consistent poten- 
tial. The access to individual nano-objects, recently de- 
velopped by different experimental techniques, provides 
a promising way of testing our theoretical results con- 
cerning the size-dependent linewidth oscillations. 

The physical relevance of the second plasmon has been 
analyzed in terms of different decay channels: Landau 
damping and particle ionization. We have shown that 
both processes are relevant, but they do not preclude 
the existence of the resonance. The comparison of our 
semiclassical calculation with the existing numerical and 
experimental results is reasonably good, despite the var- 
ious approximations of our model. 

Our theoretical results concerning the different decay 
mechanisms of the collective excitations of metallic clus- 
ters should be important for the analysis of the electron 
dynamics following short and strong laser excitations. 
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APPENDIX A: TRANSITION POTENTIAL 

In this Appendix we present the derivation of the tran- 
sition potential induced by the plasmon field and gene- 
ralize the derivation of Ref. considering the Coulomb 
interaction in the case of a dielectric mismatch between 
the electrons and the surrounding matrix in which the 
nanoparticles are embedded. Assuming that at equili- 
brium the electron density is uniform within a sphere of 
radius a, n(r) — nQ{a — r) {Q being the Heaviside distri- 
bution), a rigid displacement with a magnitude Z along 
the direction changes the density at r from n(r) to 

n(r — u) = n(r) + 5n{r). 

To first order in the field u = Ze^ , we can write 

d 

Sn{r) — — u • —n(r) = Zncos66(r — a), 
or 

We have neglected the oscillations of the density in the 
inner part of the particle due to shell effects, and also the 
extension of the electronic density outside of the particle 
(spillout effect),^ Noting Vc(r,r') the Coulomb electron- 
electron interaction, the change in the self-consistent po- 
tential due to the rigid shift (transition potential) is 

SV{r) ^ J d^r'Sn{r')Vc{r,r'). (Al) 

Using the multipolar decomposition of the Coulomb in- 
teraction, one obtains^ 

5V{r) = Z^-d(r), (A2) 

with 

{z, r a, 
za^ (A3) 
— ' ^>«- 

We notice that a displacement of the electron system 
leads to a dipolar field inside the nanoparticle, and that 
its magnitude decays as 1/r^ outside the particle. 

If we now consider the case of a noble metal nanopar- 
ticle (where the d electrons are taken into account with 
the help of a dielectric constant ed) embedded in a ma- 
trix (of dielectric constant e^), the Coulomb interaction 
between electrons is given bj^ 
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^c(r,r') 



47re^ ( 



ed f - 2/ + 1 

Lm 



r'r" (^ + l)(ed-e^) 



e„i ^ 2/ + 1 



I 



r, r' < a, 

r< ^ a, r> > a, 



r, r'>a, 



where r< = min(r, r'), r> 
the spherical harmonics. Inserting this expression into 
Eq. (|A1() , we obtain the result of Eq. HA2|I with the ad- 
ditional multiplying factor 3/(ed + 2e^). 

In both cases (with and without a dielectric mis- 
match), the expression ljA2p can be written as 6V{r) = 
ZmcUj1^d{r). The only effect of the dielectric constants 
on the transition potential as compared to the free case 
is through the red-shift of the Mie frequency. 



APPENDIX B: SEMICLASSICS WITH RADIAL 
SYMMETRY 



Semiclassical expansions constitute a very useful tool 
in mesoscopic physics since they allow for an intuitive 
description of relatively complex systems. The spectral 
properties of metallic clusters^ or the conductance fluc- 
tuations in the electronic transport through quantum 
dots^ can be readily understood when the quantum ob- 
servables are expressed in terms of an appropriate ensem- 
ble of classical trajectories. 

In problems with radial symmetry, like the one we treat 
in this work, it is tempting to take advantage of the se- 
parability into radial and angular coordinates in order 
to reduce the dimensionality of the trajectories contribu- 
ting to the semiclassical expansions. However there are 
technical difficulties introduced by the singularity at the 
origin of the centrifugal potential, and this is probably 
the reason why the radial symmetry is often not fully 
exploited in semiclassical expansions. On the other hand, 
the well-known Langer modificatior^ is a prescription 
to avoid the above-mentioned difhculties and provides 
a route to the semiclassical quantization of spherically 
symmetric systems (which has been recently extended to 
higher ordergSi). 

In this Appendix we start from the Langer modifica- 
tion in order to obtain the partial (or angular momentum 
dependent) density of states (DOS) gi{e) that we need in 
our evaluation of plasmon lifetimes. As a check of con- 
sistency, we verify in a few simple examples that when 
gi{e) is summed (in a semiclassical way) over I and m, 
we recover the well-known Berry- Tabor formula for the 
total DOS,^ 



max(r, r'), and YJ™ are 1. Langer modification and partial density of states 



For a central potential V{r), the Schrodinger equation 
is separable into angular and radial parts. The wave 
function can be written as ?/'fcim(r) = [uki{r)/r]Y/^{Q), 
where u^i verifies 



d2 



V{r) 



Uki{r) = ekiUkiir), 



2mc dr^ 

(Bl) 

with the condition Uki(0) = 0. It is important to no- 
tice that the variable r is limited to positive values and 
that the centrifugal potential possesses a singularity at 
r = 0. This significant difference between Eq. ljBl|l and 
a standard one-dimensional Schrodinger equation pre- 
vents from a naive application of the Wentzel-Kramers- 
Brillouin (WKB) approximation to treat this radial pro- 
blem. The change of variables a; = Inr and Xkiix) = 
exp {x/2)uki{r) results in a standard Schrodinger equa- 
tion. Using the WKB approximation for Xki amounts to 
change the centrifugal potential in Eq. (|B1|I according to 
the Langer modificationi2Sii£ 



l{l + l) 



1 + 



1 



The resulting WKB quantization provides the exact spec- 
trum for the hydrogen atom, as well as for the three- 
dimensional isotropic harmonic oscillator. 

The same kind of considerations in two-dimensional 
systems with a circular symmetry leads to the following 
substitution in the centrifugal potentiali^^i^ 



m 



(B2) 



which yields an exact WKB spectrum for the cases of the 
isotropic harmonic oscillator as well as for the hydrogen 
atom in two dimensions. 

The semiclassical approximation provides a method to 
calculate the leading h contributions to the DOS in the 
limit of large quantum numbers, and decomposes the 
DOS into a smooth and an oscillating part. The smooth 
term is simply the Weyl contribution— and the oscillat- 
ing term is given, in the case where the periodic orbits 
(POs) are not degenerated in action, by the Gutzwiller 
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trace formulc4i as a sum over the primitive periodic or- 
bits (PPOs). 

In the case of multidimensional integrable systems, the 
POs belonging to a torus of the phase space are degener- 
ate, and the oscillating part of the DOS is given by the 
Berry- Tabor formula as a sum over rational tori,^^ In 
one-dimensional problems, or in the radial coordinate of 
a spherically symmetric case, the trajectories are not de- 
generate, and therefore the semiclassical approximation 
to the DOS at fixed angular momentum I is given by 



with 



Qiie) 



rije) 

2TTh ' 

ttH 



osc / „\ 



(B3) 



Siie) 



where Si and t; = dSi jde are the action and period refer- 
ing to the motion in the effective (Z-dependent) radial po- 
tential; Vc (y^ is the number of classical turning points 
of the PPOs against the smooth (hard) walls. 



2. Total density of states and Berry- Tabor formula 
for systems with radial symmetry 

Using the selection rules for the plasmon decay, its life- 
time can be expressed in terms of the partial DOS Qi{e) 
whose semiclassical expression is given by Eq. (jB3|) . It is 
then important to verify that the semiclassical sum over 
angular momenta (that we use throughout our calcula- 
tions), when applied to p/(e), is able to reproduce the 
total DOS. Rather than working the most general case, 
we perform our test for three particular examples: the 
disk billiard (where the calculations are particularly sim- 
ple) , the three dimensional billiard (like the one we treat 
in the text) , and the isotropic spherical harmonic oscilla- 
tor (where the semiclassical spectrum coincides with the 
exact one). 



a. Disk billiard 



A disk billiard is defined by its radial potential 



Vir) 




(B4) 



where a is the radius of the disk. The effective radial mo- 
tion is governed by the potential V^{r) = h^m? /2mcr^ + 
V{r), with m the z component of the angular momen- 
tum included according to Eq. l|B2p . The classical PPOs 
have Vc = I'-c = \ since there is one turning point at 
the (smooth) kinetic barrier and another at the (hard) 
wall for r = a. For a given energy e we have mmax = 
{2mceY/'^a/h = (e/eo)^^^ = ka, with = 



The action and period of the PO with energy e and 
angular momentum m are given by 



m arccos 



/ m 
\ ka 



T„i{e) 



(ka)^ — m? 



ea{kaY 

respectively. The smooth part of the DOS is 



(B5a) 
(B5b) 



2TOr 



47r V 



-4, 



with A — no? being the disk area. We have replaced the 
sum by an integral and obtained the Weyl part of the 
DOS. For the oscillating part we make use of the Poisson 
summation rule and write 



+ 00 



with the phase 



m— — oo f^l 
cr=± 



2nh 



~2 



2TTmm. 



Consistently with the semiclassical expansions, we per- 
form a stationary phase approximation. The stationary 
points are given by m = ka cos ipfm, with ipfm =^ ■Kifi/f 
and the condition f ^ 2to > 0, which yield just the clas- 
sical angular momenta of the POs labeled by the topo- 
logical indices [f^fh). We then recover for the oscillating 
DOS the well-known resuHj^Si^ 



1 1 

£0 Jixka r-^ 



E/^ 



sm 



3/2 



'-Pf 



■ C0S$fm, 



ni—l f>2m 



where ffm = 1 if r = 2m and /,~m = 2 if f > 2m, 
$fm = kLfffi — 3r7r/2 + tt/A and Lfm = 2f a sin ipfm is 
the length of the orbit (f, m). 

We also notice that the quantization of the radial 
problem leads to the well-known Keller and Rubinow 
condition's, 



m arccos 



71171 + 



from which the Berry- Tabor formula can be readily ob- 
tained. 



b. Spherical billiard 

A spherical billiard is also defined by Eq. (|B4p . 
but in the three-dimensional case Vi^^{r) = h^{l + 
1/2)'^ /2'mcr^+V{r). For an energy e, the maximum value 
of the angular momentum is given by Imax = ka — 1/2. 
The action and the period of a trajectory with energy 
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e and angular momentum I are the same as in the 
two-dimensional case up to a change of to by Z + 1/2 
[Eqs. lUni]. The total DOS is given by 

^max ''max 

e^^)-Yl E Qi{e) = Y.^2l + \)gi{e). 

;=o m=~l 1=0 

For the smooth DOS, we find the first term of the Weyl 
expansioni24 



quantization rule yields the exact quantum spectrum of 
the harmonic oscillator: e„; — huj{2n + I + 3/2). 

We have demonstrated the usefulness of the radial de- 
composition for the semiclassical expansion of the DOS. 
Even in the case of degenerated classical periodic trajec- 
tories, one is able to find the semiclassical DOS by using 
the appropriate symmetry of the system, without requir- 
ing the action-angle quantization of Berry and Tabor. 



1 /2TOr 



47r2 V ?i2 



3/2 



Semiclassical dipole matrix element with 
spherical symmetry 



where V = Ana^ /3 is the volume of the sphere, and for 
the oscillating part 



1 k 



EO V TT 

^ (-l)'"sin(2y^f^)A/ cos$f^. 



f>2m 



with the same notations as in Appendix IB 2 al We again 
recover the Berry and Tabor semiclassical DOS to leading 
order in h^^^ as well as the quantization condition of 
Keller and Rubinowi^ 



c. Isotropic spherical harmonic oscillator 

The isotropic harmonic oscillator in three dimensions is 
a nonbilliard integrable system and therefore the Berry- 
Tabor quantization is very difficult to implement. The 
radial approach that we develop clearly overcomes this 
difficulty. The effective potential is V;'=^(r) = ti^il + 
l/2)2/2TOor2 -I- (l/2)mcUj'^r'^ , where u is the pulsation 
of the harmonic confinement. At a given e, we have 
^max — —1/2 + e/huj. The classical action is given by 
Si{e) = en/uj ~ TTh{l + 1/2) and the period is r n/oj. 
Using Eq. I|B3|I with i^c — 2r and v,- = Q (no hard wall) 
gives the DOS at fixed orbital momentum. 

For the smooth part of the DOS, the sum over I can be 
performed exactly, but to be consistent with the semiclas- 
sical approximation we have to take the limit e/fiio 3> 1: 
£>"(£) ~ /2{lTjidY. Writing the Poisson summation rule 
for the oscillating part and performing a stationary phase 
approximation, we have the condition on topological in- 
dices f — 2m and to ^ 1. Finally we obtain for the total 
DOS the trace formula 



2{hwf 



1 + 2 V(-l) 



m— 1 



(2iim^\ 
\ nioJ 



which has to be compared with the exact trace formula 
given in Ref. flj, where the prefactor is shifted by 
the quantity — l/8fiw, neghgible at the (high energy) 
semiclassical limit. One also notices that the WKB 



In this Appendix we focus on the semiclassical evalua- 
tion of the dipole matrix element for the case of a sphe- 
rically symmetric system, and extend the well-known re- 
sult which relates in the one-dimensional case the dipole 
matrix element to the Fourier components of the classical 
motion of the particle.'*'^ 

The spherical symmetry permits us to separate the 
dipole matrix element {nlm\z\n'Vm') into two parts: an 
angular part given by Eq. ^ and a radial part 



/ V,, 



IV 



TOc(e„ 



eni) 



dr w. 



il{r) — Ur. 

dr 



where the radial wave functions u„; satisfy Eq. (|Bip and 
we have used the commutation relation between the ra- 
dial momentum and the Hamiltonian. Next we restrict 
ourselves to the classical region in the effective poten- 
tial Vi'^^{r) between the two turning points (r_,r_|_) and 
use the WKB approximation to express the radial wave 
functions as 

2cos|l/;i^ dr'^2m, [e„, - V^^«'(^')] - 7^/4 1 



Unl{r) = 



TT{2TOe [eni-Vfir)] /ml) 



1/4 



We also assume that the radial potential is a smoothly 
varying function of the radial coordinate, that I ~ V 
(this is justified because the selection rules dictate that 
V = I ±1 and we are in the high energy limit) and that 
the energies involved in the dipole matrix element are 
sufficiently close to each other to satisfy 



2TrhAn 

Tl 



(B6) 



with Art = n' — n. 

With these approximations, changing the spatial coor- 
dinate r to the time i, we obtain, to leading order in h 



2 r'^^ ( t\ 

T^nn' = - / r(t) cos 27rAn- , (B7) 

Tl Jo \ nj 

where r(t) represents the classical trajectory in the ef- 
fective potential. Thus we see that, as in the one- 
dimensional case, the dipole matrix element is related 
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to the Fourier transform of the trajectory of the classical 
motion. 

As a check of consistency, we apply this semiclassical 
analysis to the hard-wall potential involved in our eva- 
luation of the surface-plasmon lifetime. This analysis is 
only possible in the limit sf 3> huJM '■ The approximation 
of Eq. (jB7|) is valid if we assume that the energy of the 
particle is close to the one of the hole. This energy diffe- 
rence is, because of the conservation of energy appearing 
in the Fermi Golden Rule ©, simply Hlom- 

At a given energy e, the periodic trajectory in the ef- 
fective potential is 



2TOp£ 



s$ < s$ -. 
2 



and has been approximately determined by integrating 
out the intermediate states i in the limit fcpa 3> 1. 
The integral over the intermediate state energy has been 
performed by introducing cutoffs in order to avoid un- 
physical divergencies due to the fact that discrete single- 
particle levels have been replaced in our model by a con- 
tinuum of states. When the remaining two-dimensional 
integral is evaluated numerically we obtain a smoothly 
increasing function of the parameter ^, with h(0) = as 
shown in Fig. El This function has the asymptotic limit 
lim^^oo h{^) = oo. We see that when the double-plasmon 
state is too high in energy, the linewidth diverges to in- 
finity and this resonance is no longer well-defined (the 
double-plasmon state has a lifetime equal to zero in this 
condition). 



Substituting this expression in Eq. (|B7p and making the 
expansion in 1/ An [proportional to 1/hujM, see Eq. (jB6(l ]. 
we obtain the leading order term 



7^' 



eh) 



\2 ' 



(B8) 



which agrees with Eq. Hll|) in the limit £p « Sh- We no- 
tice that this semiclassical dipole matrix element leads to 
the correct result for the smooth part r° of the single- 
plasmon linewidth in the limit ^ = hujM/£F — > of 
Eq. CEJ. 



APPENDIX C: FREQUENCY DEPENDENCE OF 
THE PLASMON LINEWIDTHS 



0.05 




FIG. 5: Functions g{^) (thick line), h{^) (full line), and 
g(5,C). The function q{^,() is represented as a function of 
^ for (^/^ = 1 (dashed line), 1.4 (dashed-dotted line), and 1.8 
(dotted line). 



In this Appendix, we present the frequency dependence 
of the single- and double-plasmon linewidths. In Fig. [S] 
we represent the function g (thick line) of ^ = S^m/ef in- 
volved in the expression of the single-plasmon linewidth 
[see Eq. (^], as well as in the first-order decay rate of 
the double plasmon (r2^i). The function g is plot- 
ted after its analytical expression [Eqs. (62) and (63) 
in Ref. and is a smoothly decreasing function with 

lim^^oo giO = 0- 

The function h involved in the expression of the second- 
order double-plasmon linewidth H31|) is defined by 

Jmax(l,25) Jo 




The function q of the two variables ^ and ( — W/ep 
involved in our evaluation of the ionization rate via the 
double-plasmon state, Eq. 1)35(1 . is defined as 



1+2^ 



dz 



(2z - 1 - C)V7"=^ 

max(2C,l+C) Z^{Z - £_){z - 1 - C) 
1 



(Cl) 



Since our approach is valid when Hujm 2fta;M, the 

function q is defined for ^ ^ C ^ 2^ and can be integrated 
numerically. The result is shown in Fig. [S] The function 
q is not very sensitive to the value of C/C = W/Hujm 
in the presented interval. However, it vanishes at the 
upper limit {W — 2hujM), since in this case particle states 



cannot be in the continuum and F;, 



0. 
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